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ABSTRACT 

We present the stellar kinematics in the central 2" of the luminous elliptical galaxy M87 (NGC 4486), using 
laser adaptive optics to feed the Gemini telescope integral-field spectrograph, NIFS. The velocity dispersion rises 
to 480 km s"' at 0.2". We combine these data with extensive stellar kinematics out to large radii to derive a black- 
hole mass equal to (6.6±0.4) x 10''Mo , using orbit-based axisymmetric models and including only the NIFS data 
in the central region. Including previously-reported ground-based data in the central region drops the uncertainty 
to 0.25 X 10^ Mq with no change in the best-fit mass; however, we rely on the values derived from the NIFS-only 
data in the central region in order to limit systematic differences. The best-fit model shows a significant increase 
in the tangential velocity anisotropy of stars orbiting in the central region with decreasing radius; similar to that 
seen in the centers of other core galaxies. The black-hole mass is insensitive to the inclusion of a dark halo in the 
models — the high angular-resolution provided by the adaptive optics breaks the degeneracy between black-hole 
mass and stellar mass-to-light ratio. The present black-hole mass is in excellent agreement with the Gebhardt & 
Thomas value, implying that the dark halo must be included when the kinematic influence of the black hole is 
poorly resolved. This degeneracy implies that the black-hole masses of luminous core galaxies, where this effect is 
important, may need to be re-evaluated. The present value exceeds the prediction of the black hole-dispersion and 
black hole-luminosity relations, both of which predict about 1 x 10^ Mq for M87, by close to twice the intrinsic 
scatter in the relations. The high-end of the black hole correlations may be poorly determined at present. 

Subject headings: galaxies: elliptical and lenticular, cD; galaxies: individual (M87, NGC4486); galaxies: 
kinematics and dynamics 



L INTRODUCTION 

The masses of central black holes in galaxies appear to be 
closely related to the luminosity (Dressier 1989; Kormendy 
1993; Kormendy & Richstone 1995; Magorrian et al. 1998) 
and stellar velocity dispersion (Ferrarese & Merritt 2000; Geb- 
hardt et al. 2000) of their host galaxies. These relationships, 
which are determined from local samples of galaxies, provide 
the means to assay the cosmological mass distribution function 
of massive black holes, and provide the empirical foundation 
for establishing the role of black holes in galaxy formation and 
evolution (e.g. Hopkins et al. 2008). 

At present the black-hole galaxy-property relationships are 
derived from several dozen black-hole mass determinations 
made over the last few decades (see Giiltekin et al. 2009). The 
relationships remain poorly observed at both their high-mass 
and low-mass ends. Lauer et al. (2007) show, for example, that 
the Mbh-ct andMsH-L relationships must be in conflict at high 
black-hole mass due to curvature in the Faber & Jackson (1976) 
relationship between galaxy velocity dispersion and luminosity. 
Small uncertainties in the high-mass end of the relations can 
lead to uncertainties of up to two orders of magnitude in the im- 
plied volume density of black holes with Mbh > 10^ Mq, due 
to the high-end exponential cutoff of the galaxy luminosity and 
velocity-dispersion distribution functions. Such estimates also 
depend critically on knowledge of the intrinsic scatter in the the 
relationships (Giiltekin et al. 2009). Thus, there remains a need 
to measure accurate black-hole masses in a sample of the most 



massive galaxies. 

Apart from the need to enlarge the sample of galaxies used to 
define the black-hole galaxy-property relationships, it appears 
that we may also need to test and potentially revise some black- 
hole mass measurements already made, especially in the mas- 
sive "core galaxies." Recent work shows that black-hole masses 
are subject to several systematic errors that have not been gen- 
erally incorporated in the models used for analyzing the data 
so far Some of these are discussed in Giiltekin et al. (2009) 
and include the radial variations in the mass-to-light ratio due 
to changes in stellar populations or the presence of a dark halo, 
uncertainties in the deprojection of the surface brightness, and 
triaxiality, among others. The most important of these system- 
atic effects are: 

Dark halo: Gebhardt & Thomas (2009) show that the mea- 
sured black-hole mass for M87 increases by more than a fac- 
tor of two when a dark halo is included in the models; the 
reason for the change is that the black-hole's kinematic influ- 
ence is poorly resolved in the data that they use, so that there 
is substantial covariance between the black-hole mass and stel- 
lar mass-to-light ratio. In-turn the best-fit stellar mass-to-light 
ratio, assumed independent of radius, is affected by whether or 
not a dark halo is included in the models. It is well understood 
that the mass-to-light profile for ellipticals changes with radius 
and not including that trend biases the black hole determina- 
tion. An obvious, but challenging, solution to this degeneracy 
is to obtain data at radii where the kinematics are strongly dom- 
inated by the black hole rather than the stars. 
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Incomplete orbit library: Shen & Gebhardt (2010) find an in- 
crease of two in the black-hole mass for NGC 4649 when using 
a more complete orbital sampling compared to models using a 
less coverage (Gebhardt et al. 2003). They argue that the or- 
bital structure near the black hole is dominated by tangential 
orbits and that the older models did not have adequate coverage 
of these tangential orbits (as discussed in Thomas et al. 2004). 
Having too few tangential orbits (i.e., too many radial orbits) 
can be compensated by having a smaller black-hole mass. 

Triaxiality: Van den Bosch & de Zeeuw (2009) find an in- 
crease of two in the measured black-hole mass for NGC 3379 
by using triaxial models compared to triaxial models (although 
they find the same black-hole mass for M32). 

All three of these systematic effects tend to increase the 
black-hole mass. The increases are generally larger than the 
statistical uncertainties and suggest that systematic effects still 
dominate. By observing stars close to the black hole, many 
model assumptions are no longer needed. For example, if the 
gravitational potential is dominated by the black hole, then 
the stellar contribution to the enclosed mass is not important; 
hence, uncertainties in the stellar mass-to-light ratio, which 
may arise from uncertainties in the dark-halo properties, can 
be mitigated by probing well inside the influence region of the 
black hole. 

Since it is among the most luminous galaxies nearby, has the 
largest black hole known (from spatially resolved kinematics) 
and has one of the nearest and best-studied AGNs, M87 is a 
natural and important target. An accurate black-hole mass de- 
termination for M87 helps to pin down the sparsely sampled 
upper end of the black-hole mass distribution and provides in- 
sights into formation and evolution of the most luminous galax- 
ies. The previous analysis of M87 from Gebhardt & Thomas is 
based on ground-based kinematic data taken in natural seeing 
under moderately good conditions (FWHM=1"). In this paper, 
we present kinematics based on the integral field spectrograph, 
NIPS, on the Gemini Telescope, taken with adaptive optics cor- 
rection. The spatial FWHM of the kinematics is 0.1" on av- 
erage, with the best seeing image at 0.08". At larger radii we 
incorporate new kinematic data out to 245" or 2.5 effective radii 
that will appear in a companion paper. The extreme improve- 
ment in the data quality of M87 allows us to model black-hole 
mass with smaller systematic uncertainty. This paper focuses 
on the determination of the mass of the central black hole; the 
analysis of the stellar mass-to-light ratio and the dark halo prop- 
erties will be given in Murphy et al. (201 1). 

Obtaining the kinematics at spatial resolution down to 0.1", 
at the same signal-to-noise obtained here, would have required 
about 100 orbits (90 hours) of Hubble Space Telescope, due to 
the faint stellar surface brightness. This adaptive optics study 
using Gemini/NIFS took about 10 hours in total, highlighting 
one of the great advantages for ground-based adaptive optics. 

We assume a distance to M87 of 17.9 Mpc. The value of the 
black-hole mass scales linearly with assumed distance. 

2. DATA 

A large amount of data exist for M87, and we do not at- 
tempt to integrate all of it; we rely on those data that provide 
the highest spatial resolution, most complete spatial coverage, 
and highest signal-to-noise. Gebhardt & Thomas (2009) com- 
bine the Hubble Space Telescope images of Lauer et al. (1992) 
with the ground-based data of Kormendy et al. (2009) at larger 
radii. These data determine the surface brightness and elUptic- 



ity from radii of 0.02" to 2000" (1.7 pc to 170 kpc). Gebhardt 
& Thomas deproject the surface brightness to obtain the stellar 
luminosity density, and we use their deprojected density pro- 
file in this paper The deprojection assumes axisymmetry; we 
assume an edge-on projection and the deprojected ellipticity is 
generally close to zero but rises in the central region to 0.2 and 
the outer region to 0.5. The stellar light profile within 0.05" 
has large uncertainties both in the radial shape and the elliptic- 
ity. The best-fit profile is a power law with exponent -0.26 in 
radius and an increase in the ellipticity inside of 0. 15". We have 
run a variety of models including and excluding this ellipticity 
change, increasing and decreasing the stellar power law with a 
range of exponents from to -0.5. We find that the black-hole 
mass changes by less than the statistical uncertainties. Thus, the 
exact shape of the central stellar light profile does not appear to 
be important for the black-hole mass. 

For the spectral data, we present new observations from the 
Gemini Telescope taken with laser adaptive optics correction 
with the integral field unit of NIFS. It is important to include 
kinematic data out to much larger radii in order to constrain 
the orbital structure and the dark halo. The main source of 
the stellar kinematics at larger radii is Murphy et al. (2011). 
They obtain spectra with the integral field unit on the McDon- 
ald Observatory 2.7m telescope (VIRUS-P, Hill et al. 2009). 
These data extend to 245". VlRUS-P has 4.1" fibers, making 
it unable to provide high spatial resolution. Thus, the region 
between the edge of the NIFS field (radius of 1.9") and 10", 
where VIRUS-P provides adequate spatially resolved kinemat- 
ics information, requires additional coverage. We therefore 
include kinematics from the SAURON integral field unit (see 
http://www.strw.leidenuniv.nl/sauron). Emsellem et al. (2004) 
present the SAURON data in terms of Gauss-Hermite polyno- 
mials. The dynamical models discussed below rely on fits to 
the line-of-sight velocity distribution (LOSVD). To convert the 
Gauss-Hermite polynomials into LOSVDs, we generate 1000 
Monte Carlo realizations of the polynomials based on their re- 
ported uncertainties. From these realizations we generate the 
average LOSVD and 68% uncertainty at each velocity bin used 
in the LOSVD. 

We only use the SAURON kinematics in this region (2.5- 
11") even though they extend to 25". We do not use SAURON 
data within 2.5" because we want to provide an independent 
measure of the black-hole mass from the NIFS data alone. In- 
cluding the SAURON data within this region does not change 
the best-fit black-hole mass, but does make the uncertainty 
smaller. We discuss these points further in §6.1. At radii be- 
yond 11", we find a difference in the dispersion between the 
SAURON values and those from the VIRUS-P data. Murphy 
et al. (201 1) argue that this difference is due to template issues 
in the kinematic extractions, and could be related to the limited 
wavelength coverage of SAURON, especially given the high al- 
pha enhancement of M87; see Murphy et al. 201 1 for a detailed 
discussion and analysis. The SAURON data used in this paper 
are available at the SAURON website. 

One option to include the SAURON data at large radii would 
be to scale the velocity dispersions in order to make the overlap 
region consistent. We do not advocate this scaling. Primarily, 
the dynamical models use the full velocity profile and not just 
moments, and it is not clear whether a simple scaling of the dis- 
persion is adequate. Furthermore, since the offset is likely due 
to template mismatch in the kinematic extraction (or continuum 
placement), the scaling may not be constant with radius. Radial 
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variations in the scaling could be due to template mix variations 
with radius, velocity dispersion changes with radius, or con- 
tinuum differences with radius due to the AGN contribution. 
Within 11" Murphy et al. find consistency with the SAURON 
kinematics. Thus, we prefer to use the SAURON data where it 
is consistent and exclude it where there are differences. 

We have performed several tests in which we exclude sub- 
sets of the data — ^removing the SAURON data, removing some 
of the large radii data, removing some of the central NIFS 
kinematics — and find no effect on the best-fit black-hole mass. 

2. 1 . Gemini NIFS observations 

We observed the central region of M87 in queue mode using 
the Near-infrared Integral Field Spectrograph, NIFS (McGre- 
gor et al. 2002), on the Gemini Telescope. We used adaptive 
optics (AO) corrections with the laser guide star system, AL- 
TAIR (Boccas et al. 2006). The AGN in M87 provided the 
low-order corrections for tip-tilt and focus in a manner similar 
to Krajnovic et al. (2009). An important feature in M87's center 
is the nearly point- like knot, located about I" off nucleus along 
the outbound jet, named HST-1 (Harris et al. 2003, Perhnan 
et al. 2003, Cheung et al. 2007, Madrid 2009), which allows 
us to monitor the telescope's point spread function (PSF). The 
PSF serves as an important input to the kinematic modeling. 
The data were taken over 5 nights in April and May of 2008 
with 23 dithered positions of 10 minute exposure each on M87. 
The telluric standards HIP 59174 and HIP 61 138 were observed 
to monitor and correct for atmospheric absorption. We used 
the K_G5605 grating that provides wavelength coverage from 
2.00-2.43 fxm, with a spectral resolution of 5290 over a field 
of 3" X 3" sampled at 0.04" north-south and 0.103" east-west 
across the image slices. 

We used the Gemini NIFS package of IRAF (Tody 1993) 
scripts (developed mainly by T. Beck) for the majority of the 
reductions. This package provides the flattening, registration 
(spatially and spectrally), hot pixel identification, and sky sub- 
traction. This package does not yet handle the error frames ap- 
propriately, so we also passed our original uncertainty frames 
(dark current, read noise, and Poisson noise) through the same 
reductions as the science data as a first modification. Our sec- 
ond modification involves interpolating over the telluric Brack- 
ett 7 line and dividing the telluric standard by a 10"* K Planck 
function; this modification to the telluric correction retains a 
proper relative spectral shape in the science data. We also 
skipped the final NIFS script step where the data are resam- 
pled to equal x and y steps since we found this interpolation 
enhanced the residual structure from our PSF fits. 

We generally took sky exposures before and after the on- 
target frames, although there were some on-target frames that 
only had one sky exposure. We used the sky nearest in time for 
each on-target exposure. The sky subtraction usually worked 
well judging from inspection of residuals under sky lines and 
tests with subtraction between different pairs of sky exposures. 
However, it is clear that some atmospheric emission line vari- 
ability occurred between our sky nods and caused uncertainty 
beyond our direct, statistical noise. There are techniques to 
bundle atmospheric emission Unes into common transition fam- 
ilies and fit a series of scalings between science and sky frames 
to minimize the residuals (Oliva & Origlia 1992, Rousselot et 
al. 2000, Davies 2007), but we have not employed them here. 
The public CO line fists extend only to 2.27 /xm and therefore 
do not cover the CO bandheads through which we measure aU 



our kinematics. However, since we have so many sets of ex- 
posures, with each set at a different dither position, the prob- 
lematic sky regions are mitigated. Furthermore, for the spectral 
extractions we mask out wavelengths near the CO bandheads 
that have large variations between sets of sky exposures above 
the thermal background. The final product from the NIFS re- 
duction package is a wavelength calibrated spectrum for each 
IFU element at each of the 23 different dither positions. We 
next find the relative position for each exposure, the PSF, and 
remove the AGN and jet continua if present. 

2.2. Determination of the PSF, Pointings and Components 

A crucial step for dynamical modeling with AO data is to de- 
termine the PSF and, in the case of M87, to remove non-stellar 
features from the spectra. For galaxies with shallow light pro- 
files like M87, determination of the PSF is particularly impor- 
tant, since one has to know how much light from the outer re- 
gions is in the central spatial elements. Fortunately, for M87 we 
are in the excellent situation of having the point source (HST-1) 
within the field that can be used as a measure of the PSF. How- 
ever, the central AGN is so bright that it significantly contami- 
nates the stellar spectra within the central few spatial elements. 

There are a few techniques to estimate the PSF with AO 
data outlined in Davies et al. (2004), which we considered. 
Davies et al. (2006) present observations where the PSF is mea- 
sured from Brackett 7 in an unresolved active galactic nucleus 
(AGN). Since M87 has Ha in the central regions, it should have 
some Bracket 7, however the redshifted line falls, to within one 
pixel, on one of the brightest sky lines at 21798A. The residu- 
als from the sky fine are strong enough to not detect Brackett 7 
emission. The only observational handles we have on the AGN 
and jet flux contributions to any particular pixel are the spa- 
tial brightness variations and the change in spectral slope. The 
AGN and jet are intrinsically much redder than the stellar pop- 
ulations. We wish to use this information without making as- 
sumptions about the stellar surface density. Thus, another goal 
is to study the stellar profile within the region dominated by the 
AGN; using the integrated fight profile does not allow one to do 
this, whereas that information is contained in the spectra. 

As opposed to measuring the light model (PSF and AGN 
fraction) from the reconstructed IFU data, we could use the 
stellar fight profile as measured from HST. There are three rea- 
sons for not forcing the light profile from a previously measured 
HST image. First, the jet has knots that are moving on short 
timescales, and we cannot be sure that the AGN and jet frac- 
tional fight and spatial position will be the same. Second, we 
desire to use the spectral information to attempt to measure the 
underlying stellar component into the center. Third, the light 
profile and AGN fraction depends on the specific color. Data 
with K band filters (Corbin et al. 2002) show color profiles that 
are flat with radius, although the analysis cannot easily go to 
radu below 1". 

We make the assumptions of a constant steUar population 
into the center, with a spatially constant spectral slope and CO 
equivalent width (EW) only altered by the relative AGN/jet 
contamination. These assumptions form a closed constraint 
on the AGN/jet components and the PSF. Similarly to Lauer 
et al. (1992) with visible light HST data and the approach of 
CLEAN algorithms (Hogbom 1974), we model the inner AGN 
jet as a set of point sources. Additionally, we fit PSFs to the 
telluric standards as verification of our in-situ PSF models. We 
find in both that the sum of an irmer, AO-corrected non-circular 



4 



Data 



HST-1 





35 40 45 50 55 



fiesidual 



200 400 800 



Fig. 1 . — Stellar/AGN/clump decomposition and PSF fit for one of the 23 datasets. In the upper left we show the original data frame collapsed across wavelength. 
In the upper right, we show the fit to the stellar distribution by enforcing a constant stellar CO EW and spectral slope across the frame. In the lower left, we show 
the fit with three point sources to the central AGN/jet, one further point source for the HST- 1 clump, and a six parameter PSF. Finally, in the lower right we show 
the residuals. The final residual map is not entirely without structure, but further point source additions do not improve the x'- Notice in the scales that different 
frames are displayed with log and linear stretched for claiity. The scales in the upper and lower left are the same (but arbitrary) units, and the scale of the residuals 
in the bottom right are in those units. Thus, the residuals are less than 1%. In this image, the inner PSF is elongated with an axis ratio of 0.6, which is can be seen 
in the left-side panels. 



function and an outer, natural seeing function fit the data well 
without spatially coherent residuals. The two-component PSF 
is common with this type of data, but circular PSFs are com- 
monly assumed (Neumeyer et al. 2007). 

It is important to get a robust spatial model for the AGN, jet 
and stellar light since the kinematics in this spectral region are 
sensitive to the equivalent width of the CO lines. Silge et al. 
(2005) show that a mismatch between the equivalent widths of 
the velocity templates with the data can bias the velocity disper- 
sion either high or low by up to 30%. Given that the additional 
continuum of the AGN will dilute the equivalent widths, it is 
important to use as much information as possible to constrain 
the relative contribution. 

Thus, we assume that i) the stellar population (and therefore 
color and spectral slope) do not vary with radius near the cen- 
ter, ii) we treat the AGN and jet as a set of point sources with 
unknown brightnesses with a flat continuum. We use the fol- 



lowing operational definition of the CO equivalent width 

£Wco = AA 2^ —T77^ (1) 



\=2.29tan 



a, X A"' 



where D{X) are the counts in each pixel, a is the fitted zero- 
point for the continuum (from a power law fit), a is the fitted 
power law exponent for the continuum, a, and as are the zero- 
point and spectral index for the stellar light. Note that under 
assumption (i) this equivalent width should not vary with posi- 
tion. 

We begin the analysis of each science frame by consider- 
ing all pixels outside of 0.9" from the center of the AGN and 
0.3" from the center of HST-1; in this way the fit to the stel- 
lar model use relatively pure stellar signal. We make a least- 
squares power-law fit from 2.1-2.27 //m to each pixel and find 
a robust estimate for EWco and as from a biweight calculation 
(Beers et al. 1990). The estimate of any pixel's stellar contin- 
uum strength, a^, can then be found with a least-squares power- 
law minimization for a,ot and a,ot followed by the application 
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of Equation 1 . The EWco and a., over all frames are 220 ± 26A 
and -3.11 ±0.43. We make continuum and equivalent width 
maps for all pixels, subtract off the non-steUar continuum, and 
normalize by the stellar continuum extrapolated through the CO 
bandheads. This procedure produces the reduced spectra which 
are later used for extraction of the kinematics. 

2.2.1. PSF Model 

The PSF determination requires further analysis. We smooth 
the stellar continuum intensities using a 0.2" x 0.2" boxcar. We 
assume that the PSF has the form of an anisotropic Gaussian 
plus a power law (in the form of a Moffatt function), given by 

PSF{x,y) = N{x,y)+M{x,y); 



Nix,y) = 



2T:a\a4 



exp(-(x^+:^)/2a^); 



M{x,y) = (1 -fliXae- l)/7ra5(l + ^4^^; 



Xc =xcos(fl3)+>'sin(a3); 
>'c=>'cos(a3)-xsin(a3), 

(2) 

where PSF{x,y) is the value of the PSF at a given position x 
and y. Nix,y) is the Gaussian model for the inner PSF, with 
normalization ai, width a2, position angle and axis ratio a4. 
M(x,y) is the Moffatt model for the outer PSF, with normaliza- 
tion (1 -ai), width 05 and exponent ag. In this model, the PSF 
is assumed constant over the NIFS 2" field; this approximation 
is very accurate for a laser guide star system at a wavelength of 
2 /xm. 

Due to the large number of parameters needed to solve for 
the PSF (shape parameters for the PSF, multiple components 
for the AGN/jet and stellar profile parameters), these fitting pro- 
cedures involve minimizing a complicated function with local 
minima. We resort to simulated annealing as a global minimiza- 
tion tool (Press et al. 1992) with temperature schedules that re- 
duce by 30% over each iteration and terminate with 10""* frac- 
tional convergence. We first perform a simultaneous fit to the 
stellar-continuum subtracted data with three central AGN/jet 
point source components, a fourth point source component for 
the HST-1 clump, and PSF parameters given by Equation 2. 

We subsample each PSF evaluation over each pixel by 5 E- 
W and 3 N-S since the individual IFU elements do not properly 
Nyquist sample the PSF. This fit determines reasonable loca- 
tions and strengths for all source components, but it improp- 
erly lets the central AGN/jet drive the PSF fit, and we know, 
in fact, that this feature is resolved given the multiple compo- 
nents. We therefore re-minimize the PSF terms and the HST-1 
terms with data within 1.2" of the preliminary HST-1 position 
but outside of 0.5" of the AGN center. The isolated clump, 
HST-1, then delivers a clean PSF determination. Finally, we 
reminimize all source components but hold the PSF terms fixed 
across the whole map. A representative decomposition exam- 
ple is shown in Figure 1 for one of the 23 datasets, where we 
show the raw IFU data, the stellar model, AGN/jet model, and 
residuals. The bottom-right panel shows the residuals plotted 
on a linear scale. There remains some structure in the residu- 
als, but the maximum residual is less than 1 % of the measured 
value. We have tried a variety of PSF models and additional 
point source models, and find no improvement. Given the tight- 
ness in the HST-1 position determination, we register all frames 



off this cleanly isolated feature. From the median of all frames, 
we estimate AGN and HST-1 spectral indices of -0.67 ±0.52 
and -1 .9 ± 2.6 respectively. 

This computation also delivers the range of PSFs for the 23 
datasets. Most of the PSFs are close to the diffraction limit of 
0.06" for the inner component, and the full range of the FWHM 
for the iimer component is 0.055 to 0.19. The fraction of light 
in the inner component ranges from 0.14 to 0.45.. The frac- 
tion of light in the central component is indicative of the strehl 
ratio; however, in practice, the strehl ratio is hard to measure 
given uncertainties in the PSF model (Gebhardt et al. 2000b). 
We make a two-dimensional image of each of the 23 analytic 
individual PSFs. From these 23 images, we then average to 
make the two dimensional array which we use for the PSF in 
the dynamical modeling. 

Figure 2 plots the flux, ellipticity and position angle versus 
radius for the combined PSF. The values are reported in Table 1 . 
This PSF has the inner Gaussian FWHM of around 0.08" with 
a fraction of the light in the central component of ai =0.38. We 
use the analysis of the PSF as measured from the reconstructed 
IFU data, and we have also compared this PSF to that as mea- 
sured from the telluric standards. We find a similar FWHM for 
the inner and outer components, but the fraction of light in the 
inner component changes. For the PSF measured from the tel- 
lurics the amount of light in the central component ranges from 
60-70%, which is 1.5 to a factor of 2 larger than that deter- 
mined from the science frames. We argue that using the telluric 
PSF is too optimistic for multiple reasons. First, the telluric 
star is used as the reference star for the PSF and tip/tilt cor- 
rections (natural guide star mode), whereas the M87 data use 
the laser as the reference. Second, the M87 data use the nu- 
cleus for the tip/tilt corrections, which is more extended than 
the telluric star. Third, the M87 frames come from long expo- 
sures of 10 minutes compare to exposure times of only a few 
seconds for the telluric. Thus, the M87 PSF is expected to be 
more extended. In any case, we also run a full set of dynamical 
models using the telluric PSF and find insignificant differences. 
It is encouraging to see that both PSFs give similar results; we 
attribute this robustness to the well-resolved kinematics in the 
region influenced by the black hole. 
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Fig. 2. — The PSF used in the dynamical modeUng. The top panel is the flux 
versus radius, the middle panel is the ellipticity profile, and the bottom panel 
is the position angle (measured N to E). The combined PSF comes from the 
sum of the 23 individual two-dimensional PSF. These values are reported in 
Table 1. 

2.2.2. Central Stellar Distribution and Offset With AGN 

Bagnuolo & Chambers (1987) and Lauer et al. (1992) both 
find power law profiles with exponent -0.26 that stays constant 
into the smallest radius measured. Subject to the constant EW 
and spectral slope assumptions of our decomposition model, we 
investigate our stellar profile and that as determined by Lauer et 
al. (1992). We find a similar slope of -0.2, and find no evidence 
for a change in the power law profile near the black hole. 

The multiple components for the AGN and jet included in 
the model provide a measure as to whether the stellar centroid 
is consistent with the AGN. Batcheldor et al. (2010) report an 
offset of 6.8 ± 0.8 parsecs using ACS images on HST. The dis- 
placement they report is along the jet axis, which they attribute 
to either a recoil event from the black hole or a black-hole 



binary. The adaptive optics data presented here have similar 
spatial resolution (0.075" for the HST data they used versus 
0.08" for the AO data). The very large advantage of the AO 
data, however, is that the spectral information provides a fur- 
ther constraint on the relative amount of AGN and stellar con- 
tribution. The average difference between the stellar centroid 
and the brightest AGN component is 3 .2 ± 6.0 parsecs, consis- 
tent with no offset. Our statistical uncertainty is 8x larger than 
that of Batcheldor et al.: 0.08" accuracy versus 0.01", respec- 
tively. Since the AO data should provide a better measure of 
the AGN and stellar contribution, it is not understood why the 
Batcheldor et al. uncertainty is 8x lower. We suspect that im- 
portant details such as the jet having multiple components (that 
may move with time) and having less observational constraints 
(imaging only versus imaging plus spectra) led to the result and 
small uncertainty reported in Batcheldor et al. We find no evi- 
dence that the AGN is offset from the galaxy center. 

2.3. Kinematic Fits 

To align the data with the kinematic axis at large radius from 
the fit of Kormendy et al. (2009), we take a position angle of 
-25° (E of N) for the major axis. Figure 3 shows the radial and 
angular bins used in the modeling, following the same birming 
of Gebhardt & Thomas (2009). Figure 3 only plots the spatial 
region for the NIFS data (the model goes out to 2000"), with a 
grayscale image from one of the 23 reconstructed IFU images. 
Spectra on opposite sides of the major axis are combined. We 
generally have ten spectra at a given radius, since we have 5 
angular bins on each side of the minor axis. The central two 
radial bins are not used due to AGN contamination (discussed 
below), and the next two outer radial bins require a sum over 
all angles in order to obtain adequate S/N. The total number of 
spectra from the NIFS data that are used for dynamical model- 
ing is then 40. Table 2 provides the locations for these 40 bins. 
We also include the signal-to-noise per pixel for each bin. 




Fig. 3 . — The binning scheme in M87 for the NIFS data only. Although this 
particular frame does not have data for each bin, the dithered set fills all bins. 
Data in the mirror bins around the major axis are added to the bins shown. 

Figure 4 show the spectra for three different spatial regions. 
The top spectrum comes from radii 0.08 <R < 0.18", where 
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we have summed over all angular bins; the middle is from 
0.18 </? < 0.3", and the bottom from a radius of 0.6". The 
wavelength range shown here is the region used for the kine- 
matic extractions. There are no significant absorption lines on 
either side, although the blueward region is used to help esti- 
mate the relative contribution of AGN and stars (as discussed 
previously). The spectrum in middle panel of Fig 4 represents 
the innermost point used in the dynamical models. The black 
line shows the data, and the dashed black lines are regions not 
used in the kinematic extractions due to large variations in the 
night sky. In this wavelength region, the residuals due to sky 
subtraction tend to be positive, even though the subtracted sky 
frames use the same exposure time as the on-target frames. We 
have tried different sky subtraction levels and find little differ- 
ence in the kinematic extractions because we excise the regions 
with sky lines. In the other wavelength regions, the sky residu- 
als average to zero. The red line in the plots is the fitted line-of- 
sight velocity distribution (LOSVD) convolved with the tem- 
plate. The template comes from a library of 10 stars observed 
with GNIRS, as reported in the GNIRS template library (Winge 
et al. 2009). We select stellar types from G dwarfs to M super- 
giants. We rely on the GNIRS template library as opposed to 
the NIFS template library since the GNIRS library contains a 
larger range of spectral types. The template library is an impor- 
tant consideration for this wavelength region (Silge et al. 2005). 
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Fig. 4. — Spectra at three different radii. The top is from 0.08 < R < 0.18", 
the middle from 0.18 < R < 0.3", and the bottom from R = 0.6". The black 
line is the spectrum and the smooth red line is the best-fit template convolved 
with the best-fit LOSVD. The dashed lines are those regions excluded from 
the fit due to high sky contamination. The spectrum at the top, which comes 
from the central region, is not used in the fit due to AGN contamination. The 
velocity dispersion obtained from the fits shown in red, from top to bottom, is 
480 km s"' , 480 km s"' , and 445 km s"' , and the S/N per pixel in each is 30, 
63, and 91. 

The kinematic extraction simultaneously fits a non- 
parametric LOSVD and template weights for the individual 
stars. The template composition is allowed to vary spatially. 
The technique is described in Gebhardt et al. (2000) and 
Pinkney et al. (2004). The LOSVD is defined in 15 velocity 
bins of 260 km s"'. There is a smoothing parameter applied 
to the LOSVD, but given the high S/N for most of the spectra, 
the smoothing has little effect on the extractions; thus, there is 
only a modest correlation between adjacent velocity bins. We 
use Monte Carlo simulations to determine the uncertainties in 
the LOSVD. The S/N of each spectrum determines the noise to 
use in the Monte Carlo simulations; from 1000 realizations of 
each spectrum, we generate an average LOSVD and the 68% 
uncertainty. 

The dynamical modeling uses the non-parametric LOSVD 
directly. However, it is sometimes convenient to express the 
LOSVD in a parameterized form as Gauss-Hermite moments, 
to show the radial run of the kinematics and to compare the data 
with the models. Table 2 shows the first four Gauss-Hermite 
moments for the NIFS data. Figure 5 plots the velocity dis- 
persion versus radius, where the dispersion is measured from a 
Gauss-Hermite fit to the LOSVDs. Figure 5 plots all of the data 
at each radius, and there are between 1 and 10 angular bins at 
each radii; thus, there are multiple points at nearly all radii in 
the figure. There is no rotation seen at a significant level in the 
NIFS data. 

We input 107 LOSVDs in the dynamical models. These 
LOSVDs come from 40 spatial bins from the NIFS data, 25 
from the SAURON data, and 42 from the large radial data of 
Murphy et al. (201 1). The data in Murphy et al. come from the 
IFU VIRUS-P, where we have nearly complete angular cover- 
age. The S/N of those data is very high (50-100 per resolution 
element). The solid line in Figure 5 plots the velocity disper- 
sion from the best-fit dynamical model. The model generates 
LOSVDs, and their dispersions come from Gauss-Hermite fits 
to those LOSVDs. For the dynamical model dispersions, we 
average along angles at a given radius for clarity. In Figure 5 
we plot both the NIFS and VIRUS-P dispersions, which have 
different PSFs. The model is convolved to each of the PSFs, 
and the plotted dispersions include the convolution. 

2.4. Spectra at R< 0.18" 

Data within the central 0.18" are excluded. Within R < 
0.08", the number of individual NIFS spatial pixels is only 50, 
whereas the number of spatial pixels for bins used in the models 
ranges from 250 to 3000. Given the shallow surface brightness 
profile for M87 and the low number of NIFS pixels, the signal 
from the central stars is low, and contamination from the AGN 
is high. We have tried a variety of models for the AGN, PSF, 
and stellar light profile; in all cases, we find that little infor- 
mation is contained in the central spectrum. We do not further 
discuss this spectrum. 

The spectrum coming from the radial region 
0.08"<R<0.18" has higher S/N but still low enough that the 
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kinematic extraction is highly uncertain. Figure 4 plots this 
spectrum in the top panel. It has lower stellar S/N compared 
to any spectra we use, and is further compromised by the un- 
certainty in the AGN subtraction. However, we still attempt a 
kinematic extraction. The red line in Fig. 4 is the best-fit con- 
volved template from the region one radial bin further in radius, 
from 0.18 < R < 0.3". This region has a velocity dispersion 
of 480 km s"' . There are wavelength regions, for example at 
2.31 fim< A <2.35 fim, where the model and data are offset. 
We do not attribute this offset to poor LOSVD modeling but 
instead to poor AGN and stellar light discrimination. We use 
this region as an example and we extracted LOSVDs including 
and excluding various regions. The range in velocity disper- 
sion over all tests is 350 to 620 km s"'. The uncertainties 
on the dispersion for each individual extraction, coming from 
Monte Carlo simulations, are much smaller than this range, 
indicating that we are dominated by systematics as opposed to 
noise. For these reasons, we exclude this spectrum. We note 
that our best-fit dynamical model predicts a velocity dispersion 
of 451 km s~' at this location (the solid line in Fig. 5 shows the 
predicted dispersion from the model). This value is in the mid- 
dle of our range of dispersions using the various extractions. 
For the central radial bin (R<0.08"), the model prediction is 
430 km s-^ 
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Fig . 5 . — Velocity dispersion versus radius for M87. The black points are the 
NIFS data. The red points are the VIRUS-P data from Murphy et al. (2011), 
and the blue points are from the SAURON data. The multiple points at each ra- 
dius represent the various position angles. The solid line is the best-fit model, 
convolved to the appropriate PSF. For the dynamical model, we include the 
predicted dispersion within 0.18". 

3. DYNAMICAL MODELS 

We use the orbit-based modeling algorithm described in Geb- 
hardt et al. (2000, 2003), Thomas et al. (2004,2005) and Siopis 
et al. (2009). These models are based on Schwarzschild's 
(1979) method, and similar models are presented in Richstone 
& Tremaine (1984), Rix et al. (1997), Cretton et al. (1999), and 
Valluri et al. (2004). 

The M87 models use the spherical geometric layout de- 
scribed in Murphy et al. (2011). We use 28 radial bins and 5 



angular bins for the spatial sampling, and 15 velocity bins. The 
smallest spatial bin goes from to 0.05". The velocity bins are 
260 km s"' wide. The average number of orbits per model is 
40,000. The orbital sampling follows the design in Thomas et 
al. (2004, 2005), and is the same as used in Shen & Gebhardt 
(2010). The latter paper illustrates the importance of a densely 
sampled orbital library: the mass obtained for NGC 4649, a 
galaxy similar to M87, is a factor of 2 larger than was found in 
earlier papers (e.g., Gebhardt et al. 2003). These papers did not 
include enough low-eccentricity polar orbits; for those galaxies 
that require significant tangential anisotropy in the central re- 
gions, this lack of circular orbits will bias the orbital structure 
and hence the black-hole mass. If a galaxy is dominated by 
tangential orbits in the central regions, the projected velocity 
dispersion will drop (for purely tangential orbits, the dispersion 
goes to zero at the galaxy center). Thus, if the central disper- 
sion is smaller than, for example an isotropic distribution, this 
drop can be accommodated by either a lower black-hole mass 
or a tangential bias with a higher black-hole mass. In fact, the 
dynamical model predicts a drop in the dispersion in the central 
region (solid line in Fig. 5); as discussed later this drop is likely 
due to a tangential bias in the orbital distribution. 
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Fig. 6. — versus black-hole mass. Each point represents the 2t that 
particular black-hole mass, and the line is a smoothed curve fitted to the points. 
The best-fit black-hole mass is (6.6 ±0.4) X lO' Mq. The vertical lines repre- 
sent the 68% range for the black hole mass. The stellar mass-to-light ratio and 
the dark halo parameters have been fixed to the values reported in Murphy et 
al. (2011). 

Gebhardt & Thomas (2009) find a strong correlation between 
the black-hole mass and the circular speed of the dark halo, 
both of which are antic orrelated with the stellar mass-to-light 
ratio. These correlations arise because of the limited spatial 
resolution of their data. This paper is based on data of higher 
quality, in particular with higher spatial resolution at the cen- 
ter (0. 1 " compared to 1 .0") and extending to larger radii (245" 
compared to 30" for the stellar kinematics). Using the present 
data, there are no significant correlations between the black- 
hole mass, stellar mass-to-light ratio, and dark halo parameters. 
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Therefore, we focus in this paper on the black-hole mass, de- 
ferring a discussion of the stellar mass-to-light ratio and dark 
halo parameters to Murphy et al. (2011). Following Mur- 
phy et al., we adopt mass-to-Ught ratio=9.1 in V (solar units) 
and a spherical dark halo with potential $(r) = jV^ log(r^ +^?)> 
Vc = 800 km s-\ Rc = 36 kpc. 

Figure 6 presents the versus black-hole mass. Each of 107 
LOSVDs that we use in the dynamical models has 15 velocity 
bins, sampling velocities from -1820 to 1820 km s~'. Given 
the dispersion profile, the outer two velocity bins at each end 
(i.e., 4 bins) are zero in the models and in the data, and since 
the uncertainties in the data are still significant for the large ve- 
locities, effectively they add nothing to x^- Thus, we have only 
1 1 LOSVD bins that contain signal (i.e., we could have limited 
the velocity range to ±1400 km s~^ with 1 1 bins and we would 
have the same x^)- The total number of data points in the kine- 
matic fits is therefore around 1 100. There is a small correlation 
between the adjacent velocity bins due to the smoothing used in 
the LOSVD extraction; this smoothing is set small enough and 
the velocity bins are large enough (260 km s~^) that the corre- 
lation only mildly reduces the number of degrees-of-freedom. 
The best-fit model has x^ = 848, so the reduced x^ is slightly 
less than unity, as expected given the correlation in the LOSVD 
bins. The points in Fig 6 are the values from the individual 
models, and the Une is a smoothing spUne. We find a black-hole 
mass of (6.6 ± 0.4) x lO'' Mq. 

Figure 7 plots the ratio of the radial velocity dispersion to 
the tangential dispersion for the best-fit model. The tangential 
dispersion is defined as = 0.5 * (f^ + cr^ + V^), where and 
9 are the spherical coordinates and is the streaming motion 
in the cj) direction. This ratio does not depend systematically on 
polar angle, and so Fig. 7 plots the angular average at a given 
radius. The confidence band comes from the range of mod- 
els that are within the 68% uncertainties of the mass model, 
based on the uncertainties of the four parameters: black-hole 
mass, stellar mass-to-light ratio, dark halo circular velocity and 
dark halo core radius. There is a sharp drop in this ratio in the 
center, implying a significant amount of tangential anisotropy 
(or similarly a lack of radial motion). As seen in Fig. 5, the 
predicted projected velocity dispersion falls strongly inside of 
0.2" (for orbits with no radial dispersion, the projected disper- 
sion in the central region would fall to zero) due to the stronger 
tangential anisotropy. At radii beyond about 30" the orbital 
structure is close to isotropic. The tendency toward stronger 
tangential anisotropy in the central region has been seen in pre- 
vious analyses for other galaxies (Gebhardt et al. 2003, 2007, 
Cappellari & McDermid 2005, Shapiro et al. 2006, Cappel- 
lari et al. 2007, Krajnovic et al. 2009). Theoretical models 
(Quinlan et al. 1995, Quinlan & Hernquist 1997, Milosavljevic 
& Merritt 2001) predict increased tangential anisotropy in the 
central regions due to a destruction of stars on radial orbits from 
ejection by or accretion onto the central black hole, leaving 
only those stars on tangential orbits. Additionally, binary black 
holes will result in a significantly increase tangential anisotropy 
(Milosavljevic & Merritt 2001), similar to the amount seen here 
in M87. While it appears that the large amount of tangential 
anisotropy seen here is due to a binary black hole, a proper anal- 
ysis requires a simulation tuned to the surface brightness profile 
and kinematics of M87. In particular, it is important to include a 
large range of initial conditions in the simulations for the stellar 
orbital structure in order to use the measured anisotropy pro- 
file to determine the evolutionary history. Given that there are 



now many galaxies with well-measured central orbital struc- 
tures, this analysis would be worthwhile. 
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Fig. 7. — Radial to tangential velocity dispersion versus radius. We average 
over polar angles in this plot since the variations in ar/ut between the angular 
bins are small. The average ratio (solid line) and 68% confidence band (dot- 
ted lines) come from all models that are within the 68% uncertainties for the 
four fitting parameters (black-hole mass, stellar mass-to-light ratio, dark-halo 
circular velocity, dark-halo core radius). An isotropic distribution would have 
the ratio equal to unity. 

4. MODELS WITHOUT A DARK HALO 

We also ran models without a dark halo to investigate the 
sensitivity of our results to assumptions about the halo. In these 
fits we include kinematic data out to a radius of 100", com- 
pared to 245" for the full dataset; we do not use kinematic data 
between 100 and 245" because in this region the kinematics 
are likely to be dominated by the dark halo. We find that the 
best-fit mass decreases to 6.4 x IQ^Mq, only 0.5-sigma or 2% 
smaller than the mass we obtain from models with a dark halo 
that use all the kinematic data. We conclude that the details of 
how we model the dark halo have negligible effect on the black- 
hole mass. In contrast, when using data with much lower spa- 
tial resolution (1" versus 0.08" in the present paper), Gebhardt 
& Thomas (2009) find a large change in the black-hole mass, 
around a factor of 2.5, between models with and without a dark 
halo. As suggested by them, the reason for this difference is 
that we now have high S/N kinematic data from well within the 
region influenced by the black hole, so the covariance between 
the black-hole mass and stellar mass-to-light ratio is negligible. 

5. M87 AND THE BH-(T AND BH-L RELATIONS 

The M87 black-hole mass derived here and in Gebhardt & 
Thomas (2009) is significantly larger than most of the the pre- 
vious determinations (with the notable exception of Sargent et 
al. 1978, which we discuss in the next section). It is thus in- 
teresting to re-evaluate M87's position in the correlations of 
black-hole mass versus velocity dispersion (BH-cr) and versus 
luminosity (BH-L). 
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5.1. The Effective Velocity Dispersion 

In Giiltekin et al. (2009) we assign M87 a velocity dispersion 
of (Tg = 375 km s"', which in turn comes from the analysis of 
Gebhardt et al. (2000). In Figure 5, however, we see that this 
value is reached only at r < 2", a location that is clearly within 
inwardly rising portion of the dispersion profile associated with 
the black-hole's kinematic influence; this value unlikely repre- 
sents the M87 galaxy overall. 

The velocity-dispersion parameter used in the BH-cr relation 
is the effective velocity dispersion, cTg, which in Gebhardt et 
al. (2000) is defined as = j^' I{r)VHr)dr / j^' I{r)dr, along 
the major axis, where / is the surface brightness and V is the 
projected second moment of the velocity distribution, and Rg 
is the half light radius, which for M87 is 100" as reported in 
Lauer et al. (2007) and Kormendy et al. (2009). This op- 
erational definition of ag appears to provide a good correla- 
tion with black-hole mass, but there are many different ways in 
which one can integrate the kinematics in order to provide one 
number for the galaxy. With this definition = 352 km s"', 
using the kinematics and surface brightness profile presented 
in this paper. The previous value of 375 km s~^ results from 
using the older kinematic and light profiles. It is clear that ag 
contains a substantial contribution from the light inside where 
the black hole influences the kinematics. If instead, we exclude 
radii within this region (defined as = GMbh / cr^ and equal to 
2.1" for our models) from the integral that determines ag, we 
find (7e = 324 km s"\ about 8% smaller We choose 324 km s~^ 
as our best estimate of (7^, with a range from 312 km s"^ to 
352kms-i. 

Churazov et al. (2010) show that there exists a radial "sweet 
spot" where the velocity dispersion at that radius is robustly re- 
lated to the circular velocity. By providing a dispersion value 
that is indicative of the galaxy as a whole, this estimate may 
correlate well with the black hole. Based on the kinematics 
from Murphy et al. (2011, the dispersion value of the "sweet 
spot" for M87 is 312 ± 10 km s~^ Furthermore, Cappellari et 
al. (2007) measure a value of 306 km s"' by integrating the 
two-dimensional data within a radius of 30". There are a va- 
riety of ways to represent a velocity dispersion for a galaxy, 
and until there is a physically-motivated model it is not obvi- 
ous which measure it optimal. Thus, in order to be consistent 
with uses of a-g for other galaxies and the current incarnation of 
the black-hole sigma correlation, one should use as reported 
above (324 km s~^), but other correlations should be studied. 

We note that contamination of dg by the light from stars 
within black hole's kinematic influence is likely to be less im- 
portant for most other galaxies, since M87 is both close and 
has an unusually massive black hole. At the same time, it may 
be prudent that this issue be examined for all galaxies in the 
context of refining the BH-u relation overall. 

5.2. Estimated Black-Hole Mass in M87 

Giiltekin et al. (2009) present two BH-ct relations, one 
for all galaxies, and one for elliptical galaxies, alone. Us- 
ing ag = 324!^! km s"^ gives \og(MBH) = 9.0!5]:2 in the first 
case, and logiMBu) = 9.1^o2 in second. Likewise, eval- 
uating the Giiltekin et al. BH-Z, relation with My = -22.71 
(Lauer et al. 2007) gives \og{MBH) = 9.0 ± 0.2. Both relations 
thus give Mbh = 1 x 10^ Mq, in contrast to our determina- 
tion of Mbh = 6.6 x 10^ Mq. Thus our measurement dif- 
fers from the predictions of this BH— cr and BH-L relations 
by 0.82 dex. However, the intrinsic scatter in these relations 



is estimated by Giiltekin et al. to be 0.44 (BH-cr, all galax- 
ies), 0.31 (BH-CT, early-type only), and 0.38 (BH-L, elUpti- 
cals only). Adding this scatter in quadrature gives estimates 
of log(MBH) = 9.0!!!;^^9.1!!|;^,9.0 ±0.4, respectively. Thus our 
measured value of \og(MBH) = 9.82 ± 0.03 differs from the pre- 
dictions by 1.4-2 sigma. Given the present uncertain state of 
knowledge of the high-mass ends of both the BH-ct and BH-L 
relations, we do remark on the larger significance of M87 devi- 
ation from the relations, except to say that is does highUght the 
need to improve the sample of black-hole mass determinations 
from the most massive galaxies. 

6. DISCUSSION 
6.1. M87 Specific Results 

Our best-fit black-hole mass is (6.6 ± 0.4) x lO'^M©. Sar- 
gent et al. (1979) report a black-hole mass of 6 x lO'M© (af- 
ter scaling to our assumed distance), which is within 1-sigma 
of our reported value. Their model is based on lower spatial 
resolution data (about 1.5"), assumes that the velocity distri- 
bution is isotropic, and does not include a dark halo. It is im- 
pressive that after three decades of improvement in data qual- 
ity, modeling, and understanding, there is essentially no change 
in the measured black-hole mass. Part of the reason for the 
robustness of the Sargent et al. result is that the radial influ- 
ence on the projected kinematics from the black hole extends to 
nearly 10" (see Fig. 5), so the influence of the black hole was 
clearly visible in their kinematic data. They also use isotropic 
models, whereas we run axisymmetric models with no restric- 
tions on the anisotropy. To study the effect of the assumption 
of isotropy, we fit isotropic models to the kinematic data pre- 
sented in this paper. The comparison between the projected 
dispersions of the isotropic models and the data is poor, with 
an increase in by over a factor of two. The poor fit makes 
it difficult to assign a best-fit mass and the range of equally 
poor fitting models have black-hole masses that range from 
6 — 8 X 10' Mq, consistent with the models of §3, which show 
significant tangential anisotropy (Fig. 7). Thus, in M87, the 
assumption of isotropy does not have a significant effect on the 
measurement of the black-hole mass, although isotropic mod- 
els provide a poor fit to the data. Sargent et al. also do not 
include a dark halo, which has been shown to cause the black 
hole to be underestimated. Their velocity dispersions at large 
radii are lower than ours (245 compared to 300 km s"'), which 
is most likely because their template library was incomplete 
and their spectra had lower S/N. The lower dispersion causes 
the assumed mass-to-light ratio of the stars to be lower, an er- 
ror of the opposite sign to the error caused by neglect of the 
dark halo. Thus, the impressive agreement between our value 
and that of Sargent et al. (1978) appears to be due in part to 
the competing effects of observational errors (dispersions too 
small, which makes the stellar mass-to-light ratio too low and 
the black-hole mass too large) and oversimplified models (no 
dark halo or velocity anisotropy, both of which make the black- 
hole mass too small). Another often-quoted black-hole mass 
determination from stellar kinematics comes from Magorrian et 
al. (1998) who report a value of 4.2 x lO^M© (for our distance). 
The likely reason for the difference is that they do not include a 
dark halo and thus overestimate the stellar mass-to-light ratio. 

The black-hole mass reported here is nearly the same as that 
reported in Gebhardt & Thomas (2009), within 4%. There is 
very little kinematic data in common between the two studies. 
The kinematic data in Gebhardt & Thomas come from older 
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long-slit data at spatial resolution of 1 .0" (van der Marel 1994), 
while in this paper we use two-dimensional coverage at spa- 
tial resolution of 0.08". We further use ground-based data from 
Murphy et al. (201 1) that have excellent S/N and radial extent. 
There is some data from SAURON (Emsellem et al. 2004) in 
common between the two studies, but this provides only 10% of 
the LOSVDs used in the models. Thus, the dynamical models 
from the two studies use nearly independent kinematic datasets, 
and give approximately the same answer. 

The uncertainties on the black-hole mass from these two 
studies are similar even though the data presented here are 
superior in many ways; the previous uncertainty is 0.5 x 10^ 
whereas the uncertainty with the NIFS data is 0.4 x 10^. In 
order to keep the black-hole mass measures independent, the 
models presented in this paper do not include the SAURON 
data inside of 2.5". The similarity in the black-hole mass un- 
certainty is then due primarily to the fact that the two sets of 
data have similar accuracy on the kinematics in the central 2.5". 
Combining all NIFS data, the accuracy on the velocity disper- 
sion is 0.2% (1 km s"^). Combining all SAURON data within 
2.5" provides the same accuracy. Thus, as long as one has a 
reUable PSF and no systematic differences in the kinematic ex- 
tractions, then it is expected that the uncertainty on the black- 
hole mass is similar using either dataset. We have run a subset 
of models including both the NIFS data and all SAURON data; 
in this case, the uncertainty on the black-hole mass decreases to 
0.25 X lO' (with no change in the best-fit mass). We report and 
utilize the result using only the NIFS data within 2.5" in order 
to 1) provide as independent result as realistically possible and 
2) control potential systematic differences in the kinematic ex- 
tractions. Murphy et al. (201 1) find a difference in the velocity 
dispersion of the SAURON data at large radii compared to their 
measurements, which they attribute to template issues. While 
we do not find an offset in the dispersion values in the central 
region, we desire to maintain the independence. The major dif- 
ference, however, is that there is no degeneracy with the stellar 
mass-to-light ratio using the NIFS data, whereas the degener- 
acy is very strong otherwise. Thus, the systematic uncertainty 
from the mass-to-light ratio profile is effectively removed with 
the adaptive optics data, making the result on the black-hole 
mass and orbital structure much more robust. 

For M87, the AO data has removed the systematics due to 
the mass-to-hght ratio profile but the systematics due to the ex- 
traction of the kinematics remain important. These systematics 
include continuum placement, template mismatch, and removal 
of the AGN contribution. The first two are general and the lat- 
ter is specific to M87. Getting any of these controlled to better 
than 1% of the velocity dispersion will be very difficult. 

Corrected to our distance, the black-hole masses reported 
from gas kinematics are (2.9 ± 0.8) X lO** Mq in Harms et al. 
(1994) and (3.8 ± 1.1) x lO^M© in Macchetto et al. (1997). As 
discussed in Gebhardt & Thomas (2009) the mass reported here 
is in conflict with these by about 2-sigma. Possible reasons for 
the differences are discussed in Gebhardt & Thomas, with the 
most likely reason being uncertainty in the inclination of the gas 
disk. Macchetto et al. assume a value of 5 1 degrees based on 
the gas kinematics. Harms et al. assume a value of 42 degrees 
based on the imaging of the gas emission. The reported differ- 
ence provides a measure of the systematic uncertainty in the in- 
chnation (i.e., whether the gas kinematics or the gas distribution 
are more affected by non-gravitational forces). Applying this 9 
degree difference in the incUnation changes the Macchetto et 



al. black-hole mass from (3.8± 1.1) to (5.4± 1.3) x 10'' Mq, 
which would lead to an insignificant difference of 0.6-sigma 
between our result and theirs. Of course, the analysis is more 
complicated than this simple application since one would need 
to re-model the gas kinematics with a different inclination. A 
proper treatment would be to include the gas kinematics with 
the stellar dynamical models. Our focus in this paper is on the 
stellar kinematics, and we do not attempt to merge the gas kine- 
matic analysis. 

6.2. General Implications for Black-Hole Mass Measurements 

While the kinematics obtained from the adaptive optics study 
produce effectively the same black-hole mass and its uncer- 
tainty from kinematics taken in native image quality, the ro- 
bustness of the measures is greatly strengthened. For exam- 
ple, the black-hole mass is not dependent on the assumption of 
constant mass-to-light ratio. Trying to generaUze this result to 
other galaxies with black-hole mass determinations is difficult 
since the measure of the black-hole mass depends on many as- 
pects. There are two observational extremes that we highhght 
as examples. The first is having a measure of black-hole mass 
that comes from observations that resolve well the kinematic 
influence of the black hole. In the most extreme case, high 
S/N spectra could potentially see the high velocity wings in the 
LOSVD due to the black hole (as discussed in van der Marel 
1994). The second example would be to allow poorer resolu- 
tion of the black hole but provide a very accurate measure of the 
mass- to-light profile. In this paper, we rely on the first strategy; 
Gebhardt & Thomas rely on the second. That the two strategies 
give consistent results, at least for M87, suggests that both may 
be reliable. 

Other studies have reported robust measures of the black- 
hole mass from ground-based studies that only poorly resolve 
the black hole's kinematic influence. Shapiro etal. (2006) mea- 
sure a black-hole for NGC 3379 from SAURON data that is 
consistent with that measured from HST data using both stars 
(Gebhardt et al. 2000c) and gas kinematics. Kormendy (2004) 
summarizes the history of black-hole mass measures for many 
galaxies and finds that, in general, the differences are within the 
reported uncertainties. If one has sufficient signal-to-noise and 
two-dimensional coverage (e.g., SAURON or VIRUS-P), then 
it should be possible to measure a black-hole mass robustly. 
Thus, it is not necessarily required to resolve the region influ- 
enced by the black hole. 

Being able to use data that does not well resolve the black 
hole's influence on the kinematics allows us to study black 
holes that are either distant or low mass. Both of these regimes 
are important for understanding the physical nature of the black 
hole correlations with the host galaxy. For example, McConnell 
et al. (2011) measure a black hole mass in NGC 6086, which 
is 133 Mpc distant. The kinematic influence of the black hole 
is barely resolved, and the degeneracy between the black hole 
mass and M/L profile is strong. However, as demonstrated for 
M87, as long as one properly characterizes the mass profile at 
large radii, then high signal-to-noise data can measure the black 
hole mass accurately. 

It is possible that systematic uncertainties bias the current 
crop of black-hole correlations. One obvious consequence 
could be that without accounting for the effect of systematic un- 
certainties, the measured intrinsic scatter would increase. Giil- 
tekin et al. (2009) measure scatter of 0.44 dex for the full sam- 
ple of galaxies with measured black-hole masses and 0.31 dex 
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for ellipticals. Once systematic effects are understood and in- 
cluded, the intrinsic scatter may decrease. Other consequences 
include increasing the mass density of black holes, if black-hole 
masses are all underestimated, and changing the slope or curva- 
ture of any correlation. Schulze & Gebhardt (2011) re-analyse 
the set of 12 galaxies from Gebhardt et al. (2003) including 
a dark halo. They find an increase of 50% in the black-hole 
mass, due primarily to improved dynamical modeling (more 
complete orbit sampling) and partly to including a dark halo. 
The increase correlates with black-hole mass. It is important to 
re-evaluate all black-hole mass estimates. 

The key to understanding all of these effects comes from 
high spatial resolution data. Data from Hubble Space Telescope 
(mainly from STIS) is generally regarded as providing the most 
significant results for black-hole mass studies. The small and 
stable PSF is a central aspect for the robustness of the data from 
HST. Future uses of STIS will play an important role for quan- 
tifying black-hole masses. The main obstacle for HST though 
is that it is a relatively small mirror and requires substantial ob- 
serving time. For example, in order to measure the black-hole 
mass in M87 at the same accuracy presented here would require 
nearly 100 orbits. While this amount of time could be justified 
for a small number of objects, going to a much larger sample 
using HST is difficult. Fortunately, adaptive optics observations 
are in a mature stage where they can provide much larger sam- 
ples. 
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